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■ Kinematic simulations of the induction equation are carried out for different setups suitable for the von-Karman-Sodium (VKS) dynamo 

■ experiment. Material properties of the flow driving impellers are considered by means of high conducting and high permeability disks that 

■ are present in a cylindrical volume filled with a conducting fluid. Two entirely different numerical codes are mutually validated by showing 
quantitative agreement on Ohmic decay and kinematic dynamo problems using various configurations and physical parameters. Field 
geometry and growth rates are strongly modified by the material properties of the disks even if the high permeability/high conductivity 

I material is localized within a quite thin region. In contrast the infiuence of external boundary conditions remains small. 

, Utilizing a VKS like mean fiuid fiow and high permeability disks yields a reduction of the critical magnetic Reynolds number for the 

■ onset of dynamo action of the simplest non-axisymmetric field mode. However this decrease is not sufficient to become relevant in the 
VKS experiment. Furthermore, the reduction of Rm"^ is essentially influenced by tiny changes in the flow conflguration so that the result 
is not very robust against small modiflcations of setup and properties of turbulence. 
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1 Introduction 

Magnetic fields of galaxies, stars or planets are produced by homogenous dynamo action in which a 
conducting fluid flow provides for generation and maintenance of field energy. During the past decade 
the understanding of the field generation mechanism has considerably benefitted from the examination of 
dynamo action in the laboratory. However, realization of dynamo action at least requires the magnetic 
Reynolds number Rm = UL/r] to exceed a threshold of the order of Rm'^'^'* ~ 10... 100. From the parameter 
values of liquid sodium - the best known liquid conductor - at standard laboratory conditions (i.e. T ~ 
200°C, r] = l/fiQa ~ O.lm^/s and C ~ Im, where fiQ is the vacuum permeability and a the electrical 
conductivity) it becomes immediately obvious that self excitation of magnetic fields in the laboratory 
needs typical velocity magnitudes of ?7 ~ lO m/s, which is already quite demanding. Therefore, the first 
successful dynamo experiments performed by Lowes and WilkinsonI (1963, 19681 ) utilized soft-iron material 



so that the magnetic diffusivity is reduced and the magnetic Reynolds number is (at least locally) increased. 
Although these experiments cannot be classified as hydromagnetic dynamos (no fluid flow and therefore 
no backreaction of the fleld on a fluid motion) they allowed the examination of distinct dynamical regimes 
manifested in steady, oscillating or reversing fields. It is interesting to note that these results did not 
initiate further numerical studies on induction in the presence of soft iron domains. 

A possibility to increase the effective magnetic Reynolds number in fluid flow driven dynamo experiments 
arises from the addition of tiny ferrous particles to the fluid medium leading to an uniform enlargement 
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of the relative permeability ( Frick et al. 2002 . Dobler et al. 20031 ) . To retain reasonable fluid properties 
the amount of particles added to the liquid is limited so that an upper boun d for the achi e vable fluid 



(|2003l l 



was 



permeability is given by /ir ~ 2. The main effect found in the simulations of iDobler et al. 
a reduced decay of the initial field but not a smaller threshold (essentially because of nonmonotonous 
behavior of the growth rate in dependence of Rm). 

Another type of ferromagnetic influence on dynamo action is observed in the von-Karman-sodium (VKS) 
dynamo. In the VKS experiment a turbulent flow of liquid s odium is driven by two counterotating impellers 
located at the opposite end caps of a cylindrical domain (iMonchaux et aLl 120071) . Dyna mo action is only 



obtained when the impellers are mad e of soft-iron with /ir ~ 100 ( Verhille et al\ 2010l ). Recently it has 



been shown in iGiesecke et al\ (|2010bl ) that these soft-iron impellers essentially determine the geometry 
and the growth rates of the magnetic field by locally enhancing the magnetic Reynolds number and 
by enforcing internal boundary conditions for the magnetic field at the material interfaces in terms of 
jump conditions. Furthermore, gradients of the material coefficients and f might support dynamo 
action because corresponding additional terms in the induction equation couple toroidal and poloidal 
fields which i s essential for the occurr ence of dynamo action. An example for this dynamo type has been 
presented in Busse and Wichtl ( 19921 ) where it was shown that even a straight fiow over an (infinite) 
conducting plate with sinusoidal variation of the conductivity is able to produce dynamo action. However, 
the experimental realization of this setup would require either an unachievable large magnetic Reynolds 
number or a rather large variation in the conductivity (^ factor of 100) whereas the mean value should 
not be too far away from the conductivity of the fluid. In order to obtain semi-homogenous dynamo action 
it might be more promising to replace the conductivity variation by a permeability variation because the 
relative permeability of soft-iron alloys easily attains values of several thousands. Although such dynamos 
are of little astrophysical relevance the experiments of Lowes and Wilkinson and in particular the rich 
dynamical behavior in the VKS dynamo demonstrate the usefulness of such models. 

The scope of the present work is to validate the numerical tool necessary to establish a basic under- 
standing of the influence of material properties on the induction process. Emphasis is given to the problem 
of the free decay in cylindrical geometry where two disks characterized by high conductivity /permeability 
and their thickness are inserted in the interior of a cylindrical container filled with a conducting fiuid. To 
demonstrate the reliability of our results we use two different numerical approaches and show that both 
methods give results in agreement. The study is completed by an application of a mean fiow as it occurs 
in the VKS experiment in combination with two high permeability disks. 



2 Induction equation in heterogenous domains 

Prom Faraday's Law in combination with Ohm's Law one immediately retrieves the induction equation 
that determines the temporal behavior of the magnetic fiux density B (often abbreviated with magnetic 
field): 

OB , IB, , , 

= \/x(uxB Vx— . (2.1) 

Ot jlQa /Xr 

In Eq. (j2.ip u denotes the fiow velocity, a the electric conductivity, /io the vacuum permeability and /ir the 
relative permeability. In case of spatially varying distributions of conductivity and permeability Eq. (j2.ip 
can be rewritten: 

dB 1 1 

_ = Vx(itxS)H ASH Vx(Vln/irXS) (2.2) 

ot AiOA*rO" AiO/^rO" 

^ -(Vln/ir + Vlncr) X (Vln/ir x B) H i — (Vln/ir + V\na) x (V x B) 



The terms on the right-hand-side that involve gradients of /ir and a couple the toroidal and poloidal field 
components which is known to be essential for the existence of a dynamo. The lack of symmetry between 
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the terms containing fi^ and a indicates a distinct impact of a and fi^ which is also manifested in the 
jump conditions for electric field and magnetic field that have to be fulfilled at material interfaces. At 
interfaces between materials 1 and 2 that exhibit a jump in conductivity a and/or in relative permeability 
fij- the normal component of the magnetic flux density is continuous whereas the tangential components 
exhibit a jump described by the ratio of the permeabilities. In case of conductivity discontinuities, the 
tangential components of the electric field are continuous and the normal component of the electric current 
is continuous. If there is no contribution of the flow, the continuity of the normal current leads to the 
discontinuity of the normal ele ctric field in th e ratio of the conductivities. Mathematically these jump 
conditions are given by (see e.g. iJacksonl . 



n-{Bi 


-B2) 


= 




_B2\ 


= 












n ■ Ui 


-32) 


= 


n X (£^1 


-E2) 


= 



(2.3) 



where n denotes the unit vector in the normal direction on the interface between materials 1 and 2. 
Although these transmission conditions are standard, their dynamical consequences in flows at large Rm 
are largely unknown. 



3 Numerical schemes 

Two different numerical algorithms and codes are used for the numerical solution of problems involving the 
kinematic induction equation (j2.ip . A combined finite volume/boundary element method (FV/BEM) is a 
grid based approach which provides a flexible scheme that utilizes a local discretization and intrinsically 
maintains the solenoidal character of the magnetic field. 

The second solution method is based on a Spectral/Finite Element approximation technique denoted 
SFEMaNS for Spectral/Finite Elements for Maxwell and Navier-Stokes equations. Taking advantage of 
the cylindrical symmetry of the domains, Fourier modes are used in the azimuthal direction and finite ele- 
ments are used in the meridional plane. For each Fourier mode this leads to independent two-dimensional- 
problems in the meridian plane. 



3.1 Hybrid finite volume /boundary element method 
We start with the induction equation in conservative form: 

dB 



where the electric field E is given by 



+Vx£; = o (3.1) 



E = -uxB + r]Vx— (3.2) 
/if 



and ?7 = 1/fioc is the magnetic diffusivity. For the sake of simplicity we give a short sketch for the 
treatment of inhomogeneous conductivity and permeability only in Cartesian coordinates. The scheme 
can easily be adapted to different (orthogonal) coordinate systems (e.g. cy lindrical or spherical coordinate 
system) making use of generalized coordinates ( Stone and Normanf 1992al lbl). 
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In the finite volume scheme the grid representation of the magnetic field is given by a staggered collo- 
cation of the field components that are interpreted as an approximation of the (cell-)face average: 



1 



AyAz 



B^{xi_,i,y,z)dydz 



(3.3) 



where the integration domain F corresponds to the surface of a single cell-face: Tyz = i , yj+i] x 
[^fc-i ' -^^fc+i] (s6e Fig. [T|). A comparable definition is applied for the electric field that is localized at the 




z 



* y 

— 



Figure 1. Localization of vector quantities on a grid cell ijk with the cell center located at {xi,yj, zi^). The dotted curve denotes the 

—i j-1 fe_i 

path along which the integration of B is executed for the computation of ^ ' ^ . 

center of a cell edge and which is defined as the line average (see Fig. [T]): 



— / E^{x,y,_i,zi,_i)d3 
Ax ' -111 



(3.4) 



Similar definitions hold for the components B^y ^ and B^'^'^ ^ , respectively for ^ '"'''^ ^ and 

2 'J 2 ■''^ 

Hi ^ 

The finite volume discretization of the induction equation reads 



dt 



^' 2'J+2. 



Ay 



Az 



(3.5) 



and it can easily been shown that this approach preserves the V • B constraint for all times (to machine 
accuracy) if the initial field is divergence free. 
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3.1.1 Material coefficients. In the following we only discuss the treatment of the diffusive part of the 
electric field, E = +rjV x B / because the advective contributions (oc —u x B) do not involve the 
mat erial properties and can be treated separately in the f ramework of an operator splitting scheme (see 
e.g. Ilskakov et qH I2nn4l . Ic^secke et oi.! 120081 . IZiegleil [l999l l . To obtain the computation directive for the 
electric field the magnetic field has to be integrated along a (closed path) around -^^(.y,^) at the edge of a 
grid cell (see dotted curve in Fig. [T]) . 



— If 1 

Ecc-r^ E^dA = - 



V X — ] dA 



V 



/\y/\z J 



(3.6) 



where F = AyAz is the surface surrounded by the path Tyz and rj is the average diffusivity {rj = [fiQa)~^) 
"seen" by the electric field. Unlike vectorial quantities the material coefficients are scalar quantities that are 
localized in the center of a grid cell. The consideration of spatial variations and/or jumps in conductivity 
respectively permeability is straightforward if corresponding averaging procedures for a or /if are applied 

to: 



respectively permeability is straigntlorward il corresponding averaging procedures lor 
( Haber and Ascheill200ll ). For the component Ex the discretization of Eq. ()3.6p leads 



e: 



B 




B,. 



(Air 



(3.7) 



In Eq. (13. 7p . r/j ■_! u_i represents the diffusivity that is seen by the electric field component E" ^ at 
the edge of the grid cell {ijk) and which is given by the arithmetic average of the diffusivity of the four 
adjacent cells: 



11i,j,k + 'ni,j-l,k + ^i-1 + f?i-l,j-l,fc 



(3.8) 



Similarly, fij- denotes the relative permeability that is seen by the magnetic field components {By and B^ 
at the interface between two adjacent grid cells which for example reads for the case given in Eq. (j3.7|) : 



for B 



ij-i,k_ 



for B 



id,k-\_ 



{l^r)i,j-\,k 
(7^r)j,i,fc-i 



2(/^r)i,j,A;(/^r)i J — l,fc 
{l^r)i,j,k + (/^r)i,j-l,fc 

2(/ir)i,j,fc(/^r)i,j,fc-l 
{fJ'v)i,j,k + (/^r)i,j,fc-l 



(3.9) 



For the computation of E 



■\,hk- 



and E^ ^ 



Eq. (1231) and Eq. 



have to be adjusted according 



to the localization and the involved field components. Applying the averaging rules (|3.8p and (|3.9p for 
the computation of the "diffusive" part of the electric field results in a scheme that intrinsically fulfills 
the jump conditions (|2.3p at material interfaces. The scheme is robust and simple to implement, however, 
the averaging procedure results in a artificial smoothing of parameter jumps at interfaces and in concave 
corners additional difficulties might occur caused by ambiguous expressions for /Xf. Furthermore in the 
simple realization as presented above, the parameter range is restricted. For larger jumps of or cr a more 
careful treatment of the discontinuities at the material interfaces is necessary which would require a more 
elaborate field reconstruction that makes use of slope limiters. 



3.1.2 Boundary conditions. In numerical simulations of laboratory dynamo action insulating boundary 
conditions are often treated adopting simplified expressions like vanishing tangential fields (VTF, some- 
times also called pseudo vacuum condition). This ambiguous denomination may be understood only for 
convecting fiows at small magnetic Reynolds number where the induced fields remain small compared to 
a given field. In fact, a restriction of the boundary magnetic field to its normal component resembles an 
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artifical but numerically convenient setup where the exterior of the computational domain is characterized 
by an infinite permeability. In case of VTF boundary conditions the calculation of field growth rates for 
typical dynamo problems are overestimated. Therefore a more elaborate treatment of the field behavior on 
the boundary is recommended which is nontrivial in non-spherical coordinate systems. Insulating domains 
are characterized by a vanishing current j oc V x S = so that B can be expressed as the gradient of a 
scalar field ^ (in the case of a simply connected vacuum) which fulfills the Laplace equation: 



B 



-V<^ with A^* = 0, ^ C»(r-^) for r — )• oo. 



(3.10) 



Integrating = and adoption of Green's 2nd theorem leads to 



^(r) 



',J G{r,r 



dn 



dn 



(3.11) 



where G{r,r') = — (47r \r — r'l)"^ is the Greens function (with AG(r,r') = —6{r — r')) and 9/dn represents 
the derivative in the normal direction on the surface element dT so that = —B^ yields the normal 
component of B on dT. From Eq. (|3.1ip the tangential components of the magnetic field on the boundary 
B^ = Br ■ B = —Br ■ V<^(r) are computed by: 



<^(r')V,^^^^ +S"(r')V,G(r,r')j dT{r') 



(3.12) 



where represents the tangential unit vector on the surface element dT{r'). After the subdivision of 
the surface V in boundary elements Tj with T = UTj the potential $i = <P{ri) and the tangential field 



BJ = BUn 



-Bj- ■ (V'Pi) in discretized form are given by 



Bf. 



(3.13) 



The solution of the system of equations p.l3p gives a linear, non local relation for the tangential field 
components on the boundary in terms of the normal components a nd closes the problem of magnetic 
induction in finite (connected) domains with insu lating boundar i es dlskakov and Dorinv 20051 ). A more 
detailed description of the scheme can be found in Giesecke et al. ( 20081 ) . 



3.2 Spectral/Finite Elements for Maxwell equations 

The conducting part of the computational domain is denoted i^c, the non-conducting part (vacuum) is 
denoted r2„, and we set Q := Qc^^v We use the subscript c for the conducting part and v for the vacuum. 
We assume that is partitioned into subregions Oci, ■ ■ ■ , ^^cAf) so that the magnetic permeability in each 
subreg ion r^c2 5 say z^*^^, is smooth. Wc denote the interface between all tlie conducting subregions. We 
denote S the interface between J7c and fly. 
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The electric and magnetic fields in 0,c and solve the following system: 



dt 

Vfi^'H'' = 0, 



V X E"", 



-u X i^'^H'' + -VxH", 
a 



di 

V-fi^'H" = 0, 
VxH'' = 0. 



V X E"-' 



(3.14) 
(3.15) 
(3.16) 



and the following transmission conditions hold across S„ and S: 



H"" xn'' + i?"' X n'' = 0, 
E;'^ X + X n'' = 0, 



(3.17) 
(3.18) 
(3.19) 



where (resp. n") is the outward unit normal on dQc H S (resp. 80,^, n S), and n^^ is the unit outward 



normal on d^ci H S^. 



3.2.1 Weak formulation. The finite element solution is computed by solving a weak form of the system 
()3.14p - ()3.19p . We proceed as follows in ilci- Multiplying the induction equation in f^ci by a test-function 
6, integrating over Qci and integrating by parts gives 



dt 



dt 



b+ [ E"-\/ xb+ [ {n" X E^') ■ b 

b+ f (-U X ^jTH"^ + -V X hA - V xb+ f E^' -{bx n^^) (3.20) 
Jn„ V ^ ) Jdn„ 



We proceed slightly differently in fi^. From ()3.16p we infer that is a gradient for a simply connected 
vacuum, i.e., = Thus taking a test- function of the form Vip, where ip is a scalar potential defined 
on r2„, multiplying (|3.14p by and integrating over Qy, we obtain 



Ft 



Vip + E"" -Vip xn" + E"" -Vi^x = 
'e Jan 



(3.21) 



We henceforth assume that a := E\qq^ is a data. Since only the tangential parts of the electric field are 
involved in the surface integrals in (|3.2Up and (|3.2ip . we can use the jump conditions (|3.19p to write 



E"' -bx n"' = / {E^} ■ b X / E'' ■ V x n'' = / E^-V^x n^ 



where {E''} is defined on by {E"} = \ {E"' + E"^) . We now add ^^TM (for i = 1, . . . , A^) and (l3:2B 
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to obtain 



dt Jn^ dt Ju£iQ„ Vo- 

+ / {E''}-lbxnj+ [ E^-{bxn^ + V^ljxn'') = - [ a-V^xn", 
Js^ Je Jan 

where we have set \b x n] := {pi x n*^* + bj x n'^-') with bi := &|n^i and bj := bln^j- We finally get rid of 
E'^ by using Ohm's law in the conductor: 

/ f ^(^7^°).v^+ f ('-V.H--u.,-H") .Vx6 (3.22) 

+ [ {-VxH" -uxfi''H''}ibxnj + [ i-VxH" -uxfi^nA ■{bxn'' + V'il^xn'') = - [ a-Vipxn" 
Je^ ct Js\o- ) Jan 

This formulation is the starting point for the finite element discretization. 

3.2.2 Space discretization. As already mentioned, SFEMaNS takes advantage of the cylindrical sym- 
metry. We denote Q^'^ and Q'^f the meridian sections of 17^ and Q^ii respectively. These sections are meshed 
using quadratic triangular meshes (we assume that O^'^ and the sub-domains ^lf.f . . . il^^ have piecewise 
quadratic boundaries). We denote {J^h}h>o^ {^h^}h>o ■ ■ ■ {J^h'^}h>o the corresponding regular families of 
non-overlapping quadratic triangular meshes. For every K in the mesh we denote Tk '■ K — > K the 
quadratic transformation that maps the reference triangle K := {(r, z) G M^, < r, < z, f + £ < 1} to 
K. Given and two integers in {1,2} with ^ we first define the meridian finite elements spaces 

:= {bh € L\n,) I bh\n.^ G CO(n,,), Vi = 1, . . . ,M, b^^TK) G P^,, G \jtx^t] , 
Xf := {^Ph G C0(O„) / i,t,{TK) G G J^} , 

where P^ denotes the set of (scalar or vector valued) bivariate polynomials of total degree at most k. Then, 
using the complex notation i^ = —1, , the magnetic field and the scalar potential are approximated in the 
following spaces: 



- 



bt{r, zy^'- Vm = 0, . . . , M, b^ G X^''"' and b^ = bl^\, 

T2=-A'/ J 

1 

^)e''"'; Vm = 0, . . . , M, e and C = V'."'" 

'i=-M ) 



where M + 1 is the maximum number of complex Fourier modes. 



3.2.3 Time discretization. We approximate the time derivatives using the second-order Backward Dif- 
ference Formula (BDF2). The terms that are likely to mix Fourier modes are made explicit. Let At be the 
time step and set := nAt, n ^ 0. After proper initialization at and t^, the algorithm proceeds as 
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follows. For n > 1 we set 



H* = IH^'"" - H^'^^-K and ( 



v,n+l ._ _ (^3^^ 



v,n+l _ ^xu,n _|_ ii;,n-l^ 



and the discrete fields g and g are computed so that the following holds for all 

C {{H'^'^+\ {b, i;)) = n{b, V), (3.23) 

where the linear for TZ is defined by 

Jan Jn,_ Jt.^ Jt. 

the bilinear form C is defined by 

£((l/'='"+\,/)'''"+i),(6,^)) := / &+ / /i^'^-^^T^ VV+ / -Vxi/^'"+i.Vx6 



At 



At 



+ g((lJ""'"+\</.^'"+^),(6,^)) + / {ivx//^'"+i}-[6xnl + /" ivxi/^'"+i- (bxn^ + V^xn^) , 

and the bilinear form (7 is defined by 

g{{Hh,iph),{bh,iJh)) ■■= Pihp^ / {Hh^ixnl + Hh^2xn2) ■ {bh,ixnl + bh,2xn''2) 

+ hhp^ {Hhxn'' + V^PhXn^) ■ {bi^xn' + V^/.xn'') , 

where hp denotes the typical size of dK U or dK U S for all K in the mesh such that dK U or 
dK U S is not empty. The constant coefficients /3i and (^2 are chosen to be of order 1. The purpose of the 
bilinear form g is to penalize the tangential jumps |i3"'^'"'"'"^xn] and H'^'^'^^ xn'^ + \/ip^'"'~^^xn^ , so that 
they converge to zero when the mesh-size goes to zero. 



3.2.4 Addition of a magnetic pressure. The above time-marching algorithm is convergent on finite 
time intervals but may fail to provide a convergent solution in a steady state regime since errors may 
accumulate on the divergence of the magnetic induction. We now detail the technique which is employed 
to control the divergence of B"^ on arbitrary time intervals. 

To avoid non-convergence properties that could occur in non-smoot h domains and discontinuous material 
properties, we have designed a non standard technique inspired from Bonito et al. ( 20101 ) to control VB. 
We replace the induction equation in VLd, i = 1, . . . ,N, by the following 



di 



-V X + i^^'Wp", (-Ao)>'* = Vti"H^\ p"\9n^, = 0. 



(3.24) 



where a is a real parameter, Aq is the Laplace operator on Q^i^ and p'^^ is a new scalar unknown. A 
simple calculation shows that p'^^ = if the initial magnetic induction is solenoidal; hence, (j3.24p enforces 
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10 



V-fj.'^^H'^^ = 0. Taking q = amounts to penalizing V^'^^iJ'^* in L'^{fici), which turns out to be non - 
convergent with Lagrange finite elements when the boundary of ilci is not smooth, (see Costabel ( 199ll ) 
for details). The mathematical analysis shows that the method converges with Lagrange finite elements 
when a € (^, 1). In practice we take a = 0.7. 
We introduce new finite elements spaces to approximate the new scalar unknown p'^ 



Xf := {p/, G L\nc) I Ph G C°(f^c.), yi = 0,...,N, phiTx) G P^,, Vi^ G ^tiJ^^', P/^ = on 917,^} 



M 



yP ■- 
- 



p 



J2 P'^ir, z)e'"'^ / Vm = 1 . . . , M, p"^ G X^'^'^ and p" 



p-^ 



m=-M 



Here ip is an integer in {1,2}. The final form of the algorithm is the following: after proper initialization, 
we solve for iJ^'^+i g , (^'''"+1 G xf and p''+^ G Xf^ so that the following holds for all b £ , ip £ 

C{{H''^+\<p-'^-^'),{b,^))+V{{H'^^+\p''^+\^^^^+'),{b,q,^)^ (3.25) 
with 

N 

V{{H,p,(l)),{b,q,^)) ■.= Y^ / (^/^"b-Vp-^"lJ"-Vg + /i2"V/^"// V^^& + /i2(i-")Vp-Vg) (3.26) 

i=i 

where h denotes the typical size of a mesh element. The term /i^" V-p'^H'^'^^^ V-^'^b is a stabilization 
quantity which is added in to have discrete well-posedness of the problem irrespective of the polynomial 
degree of the approximation for p'^. The additional stabilizing bilinear form V is defined by 

'P{(j),'il)) = I V0-V^- / ipn-Vip. 
This bilinear form is meant to help ensure that Acp'"'^^^ = for all times. 



3.2.5 Taking advantage from the cylindrical symmetry for Maxwell and Navier-Stokes equa- 
tions. SFEMaNS is a f ully n onlinear code integrating coupled Maxwell and Navier-Stokes equations 
(jGuermond et al\ (|2007l . l2009l 'l). As mentioned before, any term that could mix different Fourier modes 
has been made explicit. Owing to this property, there are M + 1 independent linear systems to solve at 
each time step (M + 1 being the maximum number of complex Fourier modes). This immediately provides 
a parallelization strategy. In practice we use one processor per Fourier mode. The computation of the non- 
linear terms in right-hand side is done using a parallel Fast Fourier Transform. Note that, in the present 
paper, we use only the kinematic part of the code with an axisymmetric steady flow. Typical time step is 
At = 0.01 and typical mesh size is /i = 1/80 with refined meshes h = 1/400 on curved interfaces. 



4 Ohmic decay in heterogenous domains 



The inspection of equations (j2.3p shows that even in the absence of flow, heterogeneous domains can lead 
to non trivial Ohmic decay problems. Therefore the reliability and the application range of both numerical 
schemes are flrst examined by studying pure Ohmic decay problems in absence of flow. A cylindrical 
geometry is chosen with radius R = 1.4 and height H = 2.6 which is in accordance with setting of the 
VKS experiment. The cylinder is filled with a conductor with diffusivity rj = (/ioO")~^ = 1 and relative 
permeability = 1. Inside the domain two disks are introduced, characterized by thickness d, conductivity 
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Figure 2. Sketch of the set up. Two disks with thickness d = 0.6, 0.3, 01 (solid, dashed, dotted curve) are introduced in a cylinder with 
height H = 2.6 and radius R = 1.4. In all runs the location of the backside of each disk is fixed at z = ±1. At the outer disk edge a 
circular shape is applied with a curvature radius corresponding to half of the disk thickness. The radial extension of the disks is fixed 
and given by i?disk = 0.95. The dashed horizontal line denotes the inner boundary that separates the dynamical active region from the 
stagnant outer layer in the runs with Rm > (see Sec.[5ll. 

a and permeability /ir (see Fig. [2|). The thickness d = 0.1 is representative of the VKS impellers but the 
other d have been tested to study the scaling law with an effective permeability or an effective conductivity 
and also to estimate the impact of the numerical resolutions. 

As long as /.tr and a are axisymmetric, in a freely decaying system the axisymmetric mode (m = 0) 
can be split into decoupled poloidal {Br,Bz) and toroidal (B^) components which decay independently 
and exhibit two distinct decay rates. The components of azimuthal modes with m > 1 are coupled and 
exhibit a single eigenstate and decay rate. In the following we limit our examinations to the decay of the 
axisymmetric mode (m = 0) and the simplest non-axisymmetric mode, the (m = l)-mode (B oc cos 99). 

4.1 External boundary conditions and field pattern 

A couple of simulations have been performed utilizing vanishing tangential field (VTF) boundary condi- 
tions (sometimes also called Pseudo Vacuum) which enforce a field geometry on the outer boundary that 
resembles the behavior in case of external materials with infinite permeability. Figure [3] shows the struc- 
ture of the field geometry with the container embedded in vacuum (upper part) and for VTF boundary 
conditions (lower part). Whereas a significant impact occurs without disks the influence of the external 
boundary conditions on the field geometry remains negligible if the disk permeability or conductivity is 
large enough. A more noticeable difference between the field distribution results from the comparison of 
high permeability disks with high conductive disks. In the first case the field structure is dominated by 
two distinct annular accumulations of azimuthal magnetic field - essentially located within the disks. The 
high conducting disk case is characterized by the domination of the axial field with a slab like structure 
concentrated around the axis in which the high conducting medium is embedded. A remarkable change 
in the field structure is obtained for the thin disk case {d = 0.1, see Fig. H] &: [5|). In case of high the 
azimuthal field is dominated by two ring like structures centered on the outer part of both disks. The 
radial field is concentrated within two highly localized paths on the outer edge of the disk whereas the 
axial field has become nearly independent from z except close to the disks because the jump conditions 
inhibit the constitution of within the disks. For the high conductive disks, the differences in the field 
pattern between d = 0.6 and d = 0.1 are less significant and the torus-like structure of the poloidal field 
component is always dominating (see right panel in Fig. [5|)' For increasing permeability, the axisymmetric 
mode changes from a poloidal dominant structure to a toroidal dominant structue (see Fig. [6] for d = 0.6). 
The mode crossing occurs for ^u^^ ~ 1.5. 

4.2 Decay rates and dominating mode 

The temporal behavior of the magnetic eigenmodes follows an exponential law B oc e"'* where 7 denotes the 
growth or decay rate. Figure [7] shows the magnetic field decay rates for a thick disk {d = 0.6) and a thin disk 
{d = 0.1) against fi*^^ (left column) and against a'^^ (right column) where p,'^^ and a*^^ denote effective values 
for permeability and conductivity that are defined as /.t^^ = J fii-{r)dV and a'^^ = J a{r)dV with 
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fioa = 100 
d = 0.6 
VTF 



Figure 3. Ohmic decay. Axisymmetric eigenmodes of the magnetic field H = B (from left to right: , H^, H^); From top to 
bottom: fi^ = Mo"" = 1 (no disks), = 100, /xgo" = 100 (all with insulating boundary conditions and d = 0.6), /ir = MO"" = 1 (no disks), 
lir = 100, iiQa = 100 (all with vanishing tangential field boundary conditions and d = 0.6). 



V the volume of the cyhndrical domain. The essential properties of the field behavior can be summarized as 
follows: The presence of high permeability/conductivity material enhances axisymmetric and (m = 1) field 
modes. However, for decreasing thickness the enhancement works selectively for the axisymmetric toroidal 
field (in case of high ^i^), respectively for the poloidal axisymmetric mode (in case of high a). For the thin 
disk the decay rate of the poloidal (respectively toroidal) field component remains nearly independent of the 
permeability (respectively conductivity). Note the changeover of a dominating axisymmetric poloidal mode 
to the dominating axisymmetric toroidal mode for the high permeability disks which occurs irrespective 
of the disk thickness around /^r ~ 1.5 (see also Fig. [6|). The mode crossing does not occur for a high 






Figure 5. Ohmic decay for thin disks {d = 0.1). Left panel: /^r = 100, right panel: fiQcr = 100). The isosurfaces present the magnetic 
energy density at 25% of its maximum value. 




Figure 6. Ohmic decay. The blue transparent isosurfaces present the magnetic energy density at 25% of the maximum value and the 
red streamlines show the field structure for d = 0.6 and (from left to right): fir = 1, 2, 10, 100 (corresponding to fil.^ = 1, 1.2, 2.7, 19.5). 

conducting disk, where the (m = 1) mode and the poloidal axisymmetric mode have nearly the same 
decay rate for /ioc''^ ^1.5 and the axisymmetric poloidal mode dominates for large a. 

Small deviations between both algorithms occur in case of thin disks (d = 0.1) for the axisymmetric 
poloidal mode and for the (m = 1) mode. A couple of simulations with higher resolution in axial direction 
(marked by the blue and the green stars in the lower right panel of Fig. [7]) show that these deviations are 
most probably the result of poor resolution in case of the FV/BEM scheme because only few grid cells 
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Figure 7. Decay rates with vacuum BC against ^'^^ (left column) and against /ioc'^** (right column) for d = 0.6 (top row) and d = 0.1 
(bottom row). The solid curves show the results obtained from the hybrid FV/BEM scheme and the dashed curves denote the results 
from the SFEMaNS scheme. The stars in the lower right panel present the results of a FV/BEM run with higher resolution demonstrating 
that the FV/BEM algorithm might approach the SFEMaNS data. 



are available to resolve the vertical structure of the disk (namely 40 mesh points for SFEMaNS against 6 
mesh points for FV/BEM). A more systematic discrepancy between both algorithms becomes obvious by 
means of the behavior of the decay time r defined by the reciprocal value of the decay rate (see Fig. [8|). For 
sufficient large (respectively a"^^), r varies following a scaling law r oc cfi'^^ (respectively cx ca^^) as 
reported in tabled) For increasing ^'^^ the decay time of the (m = 0) toroidal mode slightly increases with 
decreasing d whereas the axisymmetric poloidal mode exhibits an opposite behavior. The variation of the 
decay time with a^^ for the (m = 0) components (toroidal and poloidal) is the opposite to the ones with 
varying fi*^^. These variations suggest that the decay time scaling law is not only due to the ferromagnetic 
volume of the impellers but also to the geometric constraints associated with the jump conditions ()2.3p . 

The evaluation of the discrepancies in the scaling behavior obtained by both numerical schemes remains 
difficult because this would require larger values for /i^^ and/or a^^ which is not possible without signifi- 
cantly improving the numerical schemes. In particular for the thin disk case {d = 0.1) the achievable values 
for Hj- and/or a are restricted to iJ.f^ (respectively noa'^^) ^ 5 and, it is not obvious if the available data 
already belongs to the region that follows a linear scaling. In any case the absolute values for the decay 
rates obtained by both algorithms are close, giving confidence that the results imply a sufficient accurate 
description of the magnetic field behavior in the presence of non-heterogenous materials. 



As already indicated by the marginal differences in the field pattern for both examined boundary condi- 
tions, we find no qualitative change in the behavior of the decay rates or decay times with vacuum boundary 
conditions or VTF boundary conditions (see Fig. [9]). Although for small values of fi'^^ and a'^^ the absolute 
values of the decay rates differ by 30% the scaling behavior of the decay time is nearly independent of the 



April 2, 2010 2:4 Geophysical and Astrophysical Fluid Dynamics giesccke num'varmu 



15 



toroidal field 




(m=1)-mod. 




Figure 8. Ohmic decay. Decay times against (top row) and against fioa"^ (bottom row) for throe disk thicknesses d = 0.6,0.3,0.1 
(black, blue, green). The solid curves show the results obtained from the hybrid FV/BEM scheme and the dotted curves denote the 
results from the SFEMaNS scheme. 
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Tabic 1. Scaling coefficient c for the decay time as r oc cp.*^^^ (respectively cfiQcr^^"^) for different m — and m — 1 modes as indicated (vacuum 
BC). 



external boundary conditions (see Tab. [2|). For increasing /if the influence of these boundary conditions 
is further reduced. Whereas the decay rates (for the thick disks) differ by approximately 30% for /Xr ^ 5 
there are nearly no differences in 7 for higher values of the permeabihty. This behavior is less obvious in 
case of a high conductivity disk where the poloidal axisymmetric fleld exhibits differences in the decay 
rates of 15% even at the highest available conductivity (see Fig. llOp . Note that the axisymmetric toroidal 
field behaves exactly in the same way for both kinds of boundary conditions because insulating boundary 
conditions and vanishing tangential field conditions are identical for the axisymmetric part of B^p. 
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Tabic 2. Scaling coefficient c for the decay time as r oc cfi^ for thick disks {d — 0.6) and VTF boundary conditions. 
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Figurc 9. Decay rates and decay times against (left column) and against /jqct'^^ (right column) for vanishing tangential fields 

boundary conditions, d = 0.6 The solid (dashed) curves denote the results from the FV/BEM (SFEMaNS) scheme. 




Figure 10. Comparison of boundary conditions. Decay times against fi^ (left panel) and against fiQi^ (right panel) for vacuum BC 
(solid curves) and VTF boundary conditions (dashed curves), d = 0.6. All data results from the SFEMaNS scheme. 



5 Kinematic Dynamo 



In the following, the kinematic induction equation is solved numerically with Rm > in order to examine 
if the behavior of the magnetic field obtained in the free decay is maintained when interaction with a mean 
flow is allowed. With reference to the VKS experiment we apply the so called MND-flow ()Marie et al. 
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200i) given by 



Ur{r, z) = -{tt/H) cos{2-kz/H) r(l - rf{l + 2r), 

n^(r, z) = 4er(l — r) sin (nz/H) , 

Uz{r, z) = (1 — r)(l + r — 5r^) sin {2'kz/H) , 



(5.1) 



where = 1.8 denotes the distance between both impeher disks and e describes the relation between 
tor oidal and poloidal component of the velocity (here, e = 0.7259 is chosen following previous work, 
e.g. IStefani et alWim^ . The flow magnitude is characterized by the magnetic Reynolds number which is 
defined as 



Rm = noallmaxR, 



(5.2) 



where ?7max is the maximum of the flow velocity and a denotes the fluid conductivity. Figure [TT] shows the 
structure of the velocity field where Eqs. ()5.ip are only applied in the region between the two impellers. 
The flow active region with radius R = 1 (corresponding to 20.5 cm in the experiment) is surrounded 
by a layer of stagn ant fluid with a thickness of 0.4i? (the side layer) which significantly reduces Rm'^ 
(IStefani et alWlQOd ). n the domain of the impellers a purely azimuthal velocity is assumed given by the 
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Figure 11. Structure of the prescribed axisymmetric velocity field. The color coded pattern represents the azimuthal velocity and the 
arrows show the poloidal velocity field. The black solid lines represent the shape of the impeller disk. 



azimuthal velocity of the MND flow (Eg. I5.ip at z = Behind each impeller disk a so called lid layer 

is added. Within these lid layers a purely rotating flow is assumed, modeled by a linear interpolation of the 
azimuthal velocity at the outer side of the impeller disk towards to zero at the end cap of the cylindrical 
domain. Similar to the simulations of free decay two disks are inserted into the computational domain (see 
solid black lines in Fig. [TT]) . Here we limit our examinations to disks with a height d = 0.1. Note that the 
impellers are modeled only by the permeability and/or conductivity distribution and no particular flow 
boundary conditions are enforced on the (assum e d) inte rface between impeller and fluid. This setup is 
comparable to the configuration in iGiesecke et~ai\ (l2010bl l except that the non axisymmetric permeability 
contribution representing the blade structure has been dropped in the present study. 

Figure [12] shows the growth rates for the (m = 1) mode for different magnetic Reynolds numbers. 
Compared to the free decay, we obtain a remarkable distinct behavior of the growth rate if induction from 
a mean flow is added. A high permeability disk with Rm > causes an enhancement of the (m = 1) mode 
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compared to the case //r = 1. The shift of the growth rate increases with increasing Rm resulting in a non- 
neghgible impact on the critical magnetic Reynolds number for the onset of dynamo action. In order to get 
insight in the experimental values necessary to get dynamo action, we have computed different thresholds 
for the (m = 1) mode: Rm'^ is reduced from around 76 at = 1 to Rm'^ around 55 at /i^ = 100. The 
behavior of Rm'^ indicates a saturation around Rm*^ ~ 55 for fXj- ^ 1 which still lies above the experimental 
achievable value of approx 50. With respect to the Ohmic decay (Rm = 0) the (m = l)-mode is clearly 
suppressed with increasing fij. (see green curve in Fig. [T2]) . 

For a conducting disk we obtain a reduction of the (m = 1) growth rate. In both cases the (m = 1) 
decay rate remains independent of (respectively a) for values exceeding approximately /ir ~ 20 (or 
/UqO" ~ 20). The critical magnetic Reynolds number has also been computed for a different set-up with 




Figure 12. Growth rates for the MND flow driven dynamo against fir (left panel) and against figa (right panel). Solid curves denote 
data obtained from the FV/BEM scheme, dashed curves denote the results from the SFEMaNS scheme. The green, blue, red, yellow 
colors denote the cases Rm = 0,30, 50, 70. The black stars in the left panel show the results for the SMND flow at Rm = 50 (see text) 
as reported in Tab. [3] 



the fluid restricted to the bulk region : < r < 1.4, —0.9 < z < 0.9 with VTF conditions applied on the 
frontier of this region which results in Rm*^ = 39. Note that this pseudo- vacuum set-up under-estimates 
the threshold by more than 30%. This confirms that a realistic description of the soft iron impellers is 
crucial to get correct estimates. 

The robustness of the solutions suffers from the rather delicate dependence of the field behavior on the 
details of the flow distribution, in particular from the flow in the lid layers. Beside the well known dynamo 
killing influence of the lid flow ( Stefani et al. 20061 ) this is also true for the radial flow in the vicinity of the 
inner side of the disks. A couple of simulations have been performed applying a slightly different velocity 
field where the radial component is smoothed at the transition between the bulk of the domain and the 
impeller disk (where Ur = 0). The resulting decay rates (black stars in the left panel of Fig. [12] and Tab. [3]) 
exhibit slight differences in case of /ir = 1 and a more moderate enhancement of the (m = 1) - mode for 

= 60. 



m = 1, Rm = 50 


fir = 1(FV/BEM) 


fir = 1 (SFEMaNS) 


= 60(FV/BEM) 


fir = 60(SFEMaNS) 


MND 
SMND 


-1.218 
-1.51 


-1.327 
-1.667 


-0.550 
-1.16 


-0.655 
-1.291 



Table 3. Decay rate for m — 1 mode for 2 flows MND and a similar flow with slightly modifled (smoothed) radial velocity component (SMND). 
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6 Conclusions 

Although soft iron is strongly connected to magnetostatics, experimental dynamos have shown that this 
material may also find important applications in the field of magnetohydro dynamics. For instance, at least 
one of the two impellers of the Cadarache experiment must be made of soft iron in order to achieve dynamo 
action. This is an unexplained fact which suggests that one may wonder if the role of this material is only 
to lower the critical magnetic Reynolds number in the domain of experimental feasibility or if the dynamo 
mechanism is fundamentally different when the conducting medium is no longer homogenous. This issue 
may be faced in principle numerically. However, to face such problems with heterogenous domains, specific 
algorithms must be implemented and validated and this is the aim of the present study since analytical 
results are lacking. Our comparative runs of Ohmic decay problems proved in practice to be extremely 
useful to optimize both codes and to select some numerical coefficients occurring in the algorithms (such 
as in penalty terms). 

The problems which have been successively presented above are standard in MHD, but we were forced 
to reduce the dimension of the parameter space to configurations more or less related to the Cadarache 
experiment, where the impellers may be treated as disks in a conducting flo w bounded by a cylinde r of 
a given aspect ratio. We have thus considered axisymmetric domains only (see Giesecke et al. ( 2010bl ) for 
non-axisymmetric cases), and azimuthal modes of low order (m = and 1). 

We have first studied Ohmic decay problems, with disk impellers of various thicknesses to investigate 
scaling laws and the impact of the spatial resolution. Internal assemblies of high permeability material 
within the fluid container are different from the problem of an enhanced, but homogenous fluid perme- 
ability because of inner boundary conditions for the magnetic field (in case of high permeability material), 
respectively for the electric field/current (in case of conductivity jumps). In the free decay problem with 
thin high permeability disks a selective enhancement of the axisymmetric toroidal field and the [m = 1) 
mode is observed whereas the axisymmetric poloidal field component is preferred in case of high conductive 
disks. 

We have also shown that pseudo-vacuum boundary conditions, which are easier to implement on the 
cylinder walls than the jump conditions on the impellers, have only a slight influence on the decay rates. The 
impact of the outer container boundaries on the field behavior is limited to a shift of the decay/growth 
rates. This is surprising, insofar as pseudo vacuum boundary conditions resemble the conditions that 
correspond to an external material with infinite permeability. Nevertheless, the presence of high perme- 
ability/conductivity disks within the liquid occlude the influence of outer boundary conditions, and the 
simplifying approach applying vanishing tangential field conditions at the end caps of the cylinder in order 
to mimic the effects of the high permeability disks in the VKS experiment is not sufficient to describe 
the correct field behavior. The consideration of lid layers and disks with (large but finite) permeability 
remains indispensable in order to obtain the influence of the material properties on growth rates as well 
as on the field geometry. 

For completeness, we have also considered conductivity domains. Prom the experimental point of view 
the utilization of disks with a conductivity that is 100 times larger than the conductivity of liquid sodium 
remains purely academic. Nevertheless, the simulations show a crucial difference between heterogeneous 
permeabilities and conductivities: even if these two quantities may appear in the definition of an effective 
Reynolds number Rm^''^ = fiQ^i^^ a'^^U L, they do not play the same role and they select different geometries 
of the dominant decaying mode. It is not only a change of magnetic diffusivity that matters. 

We considered then kinematic dynamo action, using analytically defined flows in accordance with the 
setting of the VKS mean flow. Since these flows are axisymmetric, the azimuthal modes are decoupled. 
The most important is the (m = 1) mode which will be excited eventually through dynamo action. We 
have shown that our codes give comparable growth rates for this mode. We have examined also the growth 
rate of the (m = 0) magnetic field in presence of soft iron impellers and the axisymmetric MND flow. Since 
convergence of results is not achieved in all the cases considered, this comparative study is still in progress 
and it has thus not been included in the present paper. We recall that the main surprise of the Cadarache 
experiment was perhaps the occurrence of the mode (m = 0), which pointed out the possible role of the 
non-axisymmetric flow fluctuations. Non-axisymmetric velocity contributions might be considered in terms 
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of an a-effect as it has been proposed in lPetrelis et all ( 20071 ) andlLague rre et~al. (l2008a b ). Preliminary 
examinations applying simple a-distributions are presented in LGlcsccke e t al\ (l2010aF and iGiesecke et ali 
(l2010bl V ^lowever, there is still a lack of knowledge on the details and physical justification on a precise 
a-distribution which requires a non-linear hydrodynamic code. The questions related to this empirical 
fact represent a main issue of the experimental and numerical approaches of the fluid dynamo problem 
and deserve a dedicated study. So far our model is not capable to explain the main features of the VKS 
experiment, which are the dominating axisymmetric field mode and the surprising low critical magnetic 
Reynolds number of Rm ~ 32. However, our results give a hint why the (m = 1) mode remains absent in 
the experiment. Dynamo action may occur when coupling terms between the magnetic field components 
are present and antidynamo theorems are derived when such terms are lacking: this is in particular the case 
for the (m = 0) mode with an axisymmetric flow (Cowling's theorem). Conversely, a source term on this 
mode appears when the flow axisymmetry is broken. Although the relative amplitude of this source cannot 
be discussed here, we note that the decay time of the (m = 0) toroidal mode become the smallest when 
the effective permeability is high enough (see for example Fig. [Tj). It may thus appears as the dominant 
mode of the dynamo, as it seems to be observed in the VKS experiment. Otherwise stated, the impact 
of soft-iron impellers on the critical magnetic Reynolds number of the (m = l)-mode could be rather low 
(decrease from 76 to 55 in the MND case) and could remain unobservable, while it could be strong 
for the (m = 0) mode (down to 32 in the VKS geometry) when conjugated to a slight departure from flow 
axisymmetry. Numerical evidences of this picture are growing. 
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